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Patterns in reaction-diffusion systems often contain two spatial scales; a long scale 
^ , determined by a typical wavelength or domain size, and a short scale pertaining to front 
^ ! structures separating different domains. Such patterns naturally develop in bistable and 



excitable systems, but may also appear far beyond Hopf and Turing bifurcations. The 
global behavior of domain patterns strongly depends on the fronts' inner structures. In 
this paper we study a symmetry breaking front bifurcation expected to occur in a wide 
class of reaction-diffusion systems, and the effects it has on pattern formation and pattern 
dynamics. We extend previous works on this type of front bifurcation and clarify the 
relations among them. We show that the appearance of front multiplicity beyond the 
bifurcation point allows the formation of persistent patterns rather than transient ones. In 
a different parameter regime, we find that the front bifurcation outlines a transition from 
oscillating (or breathing) patterns to traveling ones. Near a boundary we find that fronts 
beyond the bifurcation can reflect, while those below it either bind to the boundary or 
disappear. 
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1. Introduction 

Pattern formation out of equilibrium can often be attributed to the coupling of dif- 
fusion with local nonlinear dynamics. These processes are explicit in reaction-diffusion 
systems such as chemical reactions or electrical activation of biological membranes, and 
implicit in systems whose large scale behavior is governed by dissipative amplitude equa- 
tions such as lasers. The coupling of nonlinearity and diffusion is most beautifully realized 
in two chemical systems; the Belousov-Zhabotinsky (BZ) reaction, where a variety of 
traveling patterns have been observed (for a recent review see [Mer92]), and the Chlorite- 
lodide-Malonic-Acid (CIMA) reaction, where stationary, traveling, and turbulent patterns 
have been found [CDB90,OuS91,LKE92]. Reaction-diffusion systems can be divided into 
three main categories: Excitable, bistable, and Hopf- Turing systems. Patterns in excitable 
and bistable systems normally involve front structures separating domains of different uni- 
form or quasi-uniform states. They coexist with one or two stable uniform states and thus 
need to be triggered by specific choices of initial conditions. A Hopf- Turing system, on 
the other hand, gives rise to patterns through the destabilization of a uniform stationary 
state in a Hopf or a Turing bifurcation [Tur52,Mur89]. Close to onset, smooth oscillatory 
or stationary patterns usually appear. However, as the system is driven away from on- 
set, short scale structures, in the form of boundary layers or fronts, may develop. In this 
regime, the patterns that emerge resememble those in excitable and bistable systems. 

In this paper we study patterns in reaction-diffusion systems consisting of fronts sep- 
arating different domains. Our primary concern here is with the effect front bifurcations 
may have on the dynamics of fronts and on the patterns they form. Most of our analysis 
will focus on bistable systems, however, many of the results to be derived also apply to 
excitable systems, and to Hopf- Turing systems far beyond onset. 

We shall consider a specific reaction-diffusion model that exhibits the same variety of 
patterns that have been observed in both the BZ and CIMA reactions. It shares many 
gross features with models of the BZ reaction [Fib85] and the Lengyel-Epstein model of 
the CIMA reaction [LeE92], but is more amenable to analysis. The model consists of 
two scalar fields, u{x,t) and v{x,t) and contains four parameters; the ratio, e = T^/T^, 
between the time scales associated with the two fields, the ratio, S = D^/Du, between 
the two diffusion constants, and two parameters, ai > and ao, characterizing the local 
reaction dynamics. The model reads 

u^ = u - - V -\- Uxxi (1.1a) 



vt = e{u - aiv - ao) + Sv^x, 



(1.16) 
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where the subscripts on u and v denote partial derivatives. Models of this type have 
extensively been studied in the context of excitable media where e <^ 1. The form with 
(5 = is known as the FitzHugh Nagumo model of nerve conduction [Fit61,NAY62]. In the 
following we will refer to (1.1) as the FitzHugh Nagumo (FHN) model for (5 7^ as well. 

The stationary homogeneous states of (1.1) are determined by the intersection points 
of the two nuUclines v = u — and v = {u — ao)/ai. Three basic cases, corresponding 
to the three types of reaction-diffusion systems mentioned above, can be distinguished as 
Figs. 1 illustrate, (z) The nuUclines intersect at a single point lying on one of the outer 
branches of the cubic nuUcline v = u — (Fig. la), (n) The nuUclines intersect at a 
single point lying on the middle branch of the cubic nuUcline (Fig. lb), {iii) The nuUclines 
intersect at three points each lying on a different branch of the cubic nuUcline (Fig. Ic). 
Note that these intersection points, and consequently the stationary homogeneous states 
of (1.1), are independent of the parameters e and 5. These parameters, however, do 
affect the stability of the stationary homogeneous states, and will serve in the following 
as bifurcation parameters. They also have physical meanings that do not depend on the 
details of the FHN model. Thus, with appropriate choices of e and 5, case (z) can describe 
an excitable medium, case (ii) a medium undergoing either an Hopf or a Turing bifurcation, 
and case {iii) a bistable medium. Note also that when oq = 0, Eqs. (1.1) have an odd 
symmetry {u,v) — > {—u,—v). Analyzing the symmetric system will help us identifying 
symmetry breaking bifurcations and multiplicity of stable front solutions. Finally, a few 
words about notations. The solutions u{x,t) and v{x,t) of (1.1) obviously depend on the 
four parameters, qq, ai, e, and S, but we will avoid displaying this dependence unless we 
specifically wish to address it. In that case we will display only the parameter (s) under 
consideration (see for example (3.25)). The same rule holds for critical values of e or 5 at 
which bifurcations occur (see (3.33)). 

We begin in section 2 by first identifying Hopf and Turing bifurcations for case (ii). We 
then use singular perturbation considerations to argue that far beyond these instabilities 
front structures can develop, and that these structures coincide, to leading order, with 
the single front solutions of case (Hi). This will allow, later on, extending some of the 
conclusions we draw for bistable systems to Hopf- Turing systems driven far beyond onset. 
In section 3 we proceed to case (Hi) or bistable systems and study a symmetry breaking 
front bifurcation that gives rise to front multiplicity. This bifurcation is closely related 
to chiral symmetry breaking in magnetic domain walls [LaN79] and in front structures 
arising in forced oscillatory systems [CLH90]. Indeed, some of the results to be derived 
here apply to these systems as well. The apearance of front multiplicity has dramatic 
effects on the behavior of patterns. For 5 sufficiently small it allows for the formation of 
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persistent patterns (stable traveling waves) rather than transient ones, whereas for 6 large 
enough it is responsible for transitions from oscillating or stationary domains to traveling 
ones. Another consequence of the coexistence of multiple fronts is the possible reflection of 
traveling fronts at boundaries. We study these effects in section 4. In section 5 we describe 
the numerical procedures we have used in our simulations. We conclude in section 6 with 
a brief summary and a discussion of new problems motivated by the present work. 

2. The Development of Fronts Far Beyond Onset 

The FHN model, like a number of other reaction-diffusion models [RoM92], may 
undergo both Turing and Hopf bifurcations. Figs. 2 show the patterns that develop near 
and far from the onset of these instabilities. The appearance of patterns involving two 
spatial scales far beyond onset can be understood using singular perturbation theory as 
we now show. 

For the sake of simplicity, consider the symmetric FHN model, that is (1.1) with 
ao = 0, and assume further that < ai < 1. These parameter settings imply case (ii) 
(see Fig. lb). The stationary homogeneous state, (u, v) — (0, 0), loses stability in a Turing 
bifurcation when 

^-1 = = 2 - ai + 2(1 - ai)^/^ and e > 1/ai, (2.1) 

where 

fi = e/S. 

For jJT^ > 11^^ the stationary homogeneous state becomes unstable to perturbations of 
finite wavenumbers. The same state loses stability in a Hopf bifurcation when 

e = eH = l/ai. (2.2) 

For e < it becomes unstable to uniform perturbations. 

Far beyond the Turing bifurcation ji = e/5 <^1. Introducing a rescaled space coordi- 
nate z = ^/Jix., we find that stationary solutions of (1.1) should satisfy 

l^-Uzz -\- u — V — v?' — (2.3a) 

f^zz -\-u — aiv = 0. (2.36) 

In all regions where u varies on the same scale as v we can neglect the second derivative 
term in (2.3a) and solve the remaining cubic equation for u in terms of v. There are three 
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solution branches of which the outer two, u = U- (v) and u = u+{v), represent, respectively, 
low and high u values. In the following we will call regions of high and low u values up state 
and down state domains, respectively. Using u±{v) in (2.36) we obtain closed equations 
for V in the up and down state domains. There also might be regions where the second 
derivative term in (2.3a) balances the remaining terms. In these regions u varies on a scale 
of order ^//J which is much shorter than the scale of order unity over which v varies. Such 
boundary layers arise when nearby regions converge to different branches u±{v) (because 
of appropriate initial conditions). Stretching back the space coordinate we find to leading 
order: 

"^xx + u — V — u^ = Vxx = 0. (2-4) 

These equations have two symmetric front solutions 

v = 0, u = ±taiih{x/V2). (2.5) 

To construct approximate solutions representing periodic arrays of front structures as in 
Fig. 26, we solve the closed equations for v in the up and down state domains, Vzz + 
U'±{v) — aiv = 0, and match the solutions and their first derivatives at the front regions. 
For short wavelength patterns we may simplify the equations for v by linearizing u±{v) 
around v = 0: u±{v) ~ ±1 —v/2. Carrying out this calculation we find a band of periodic 
stationary solutions with up and down state domains of equal size, a consequence of the 
odd symmetry of the model when qq = 0. For more details about this type of calculation 
the reader is referred to Refs. [DKT88,KeL89,Mer92] as well as to Appendix B, where the 
size of a single domain structure in an asymmetric bistable system is calculated. 

Similar considerations apply to the Hopf bifurcation case. Far beyond that bifurcation 
e <S ch, or e <^ 1 if we choose ai ~ 0{1). Traveling wave solutions of (1.1) satisfy 

ew^^ + ecou^ + u — V — = 0, (2.6a) 

Sv^(^ + cqv^ + u — aiv = 0, (2.66) 

where ( :— \^e{x — ct) and c = coVe is the traveling wave speed. This choice of c is 
appropriate for short wavelength patterns [DKT88], but other choices may work as well 
[Fif85]. As before, we distinguish between domains where both u and v vary on a scale of 
order unity, and fronts where u varies on a scale much shorter than that of v, this time 
of order ^/e. For domains we find again the relations, u = u±{v), and obtain the closed 
equations for v: dv(^(^ + cqv^^ + u±{v) — aiV = 0. To study the fronts we stretch the moving 
coordinate frame, x — C/v^ — x — ct, and find the leading order equations 
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These equations have the front solution 

v^Vf, [u_(w/)e«^ + u+(i;/)]/[l + e«^], (2.8) 

where q = [u+{vf) — U-{vf)]/y/2 and Vf is the (approximately) constant value of v across 
the front. A symmetric front solution is obtained from (2.8) by the transformation x — > —x. 
Periodic traveling solutions as appear in Fig. 2d can be constructed as before. 

Note that the singular perturbation analysis sketched above does not depend on the 
parameter ai in any way that would affect the qualitative results (at least for patterns 
whose wavelengths are suflBciently short). In fact, for the particular model considered 
here, the front structures are strictly independent of ai (see (2.5,8)). We therefore expect 
to find similar patterns in bistable systems for which ai > 1. 

3. A Front Bifurcation 

Despite the vast literature on FHN type models very few studies addressed other than 
small e values. The significance of broadening the range of e is that there exists a critical 
value ec at which a front bifurcation occurs. This bifurcation has important implications 
on pattern formation as we will see in sections 4 and 5. Rinzel and Terman [RiT82], who 
studied the FHN model with 5 = 0, were the first to observe the creation of fronts at a 
critical e value. More recently, Ikeda et al. [IMN89] analyzed the regime 5/e » 1, and 
found a pitchfork front bifurcation. In this section we extend these results in a number 
of ways. First, we study the FHN model with 5 = and obtain explicit forms for the 
front solutions and for the bifurcation diagram in the vicinity of the bifurcation point. 
Then we consider nonzero 5 values and evaluate a bifurcation line, 5 = 5c{e), in the e — 5 
plane. Finally, we connect these results to a recent work by CouUet et al. [CLH90] where 
a nonequilibrium analog of the transition from Ising to Bloch walls in ferromagnets with 
weak anisotropy [LaN79] has been found. This relation gives a broader perspective to the 
body of works on patterns in FHN type models. Before embarking on the front bifurcation 
analysis we discuss a few limiting cases of (1.1) which motivate some aspects we wish to 
emphasize later on. 

3.1 The variational case: e = 

The simplest front structures appearing in dissipativc systems are those involving a 
single scalar field. One example, known in the context of equilibrium phase transitions, is 
the time dependent Ginzburg-Landau (TDGL) model obtained from (1.1) by either setting 
e = 0, or regarding f to be a constant parameter: 

Ui. = u-U^ -V + Uxx- (3.1) 
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This equation has the variational form 

ut = -ST/Su, (3.2) 

where J-' is a Lyapunov functional given by 

r = J [£{u,v) + ul/2]dx, £{u,v) ^ -u^/2 + u^/4: + vu. (3.3) 

In the terminology of equilibrium phase transitions, u is an order parameter, v is an 
external field, and .F is a free energy. For v values in the range -2/(3^3) < ^; < 2/(3x/3) 
the free energy density, has double-well forms as shown in Figs. 3. The two wells 
correspond to the stationary homogeneous states, U-{v) and w+(f), that solve the cubic 
equation — u + v = 0. 

Front solutions of (3.1) (or domain walls separating different phases) propagate in a 
preferred direction dictated by the minimization of T. The speed of propagation, c, is 
determined by the nonlinear eigenvalue problem 

u" + cu' + u-v-u^ = 0, (3.4a) 

-(X) = i"+S1 ^^""^ (3.46) 

where x = ^ ~ ct. An expression for c can be obtained by multiplying (3.4a) by u' and 
integrating along the whole line. This yields 

c — c{v) — a{v) £(u-{v)) — £(u-^{v)) , (3.5) 

where a{v) = 1/ u'{v)^dx is positive. For v = Q the two wells of £ are of equal depth 
and the front solution of (3.4) is stationary (c = 0). For negative v values, £{u-) > £{u+) 
and the front moves in the positive x direction (c > 0) so as to increase that part of the 
system having lower energy. When v is positive the front propagates toward negative x 
values (c < 0). The explicit form of this front solution is given in (2.8). The symmetric 
solution is obtained by the transformation x ~^ ~X ^^'^ c — > — c. 

3.2 The regime e > 1 

We turn now to the full form of (1.1) with 5 = 0, and study the case of bistability {in) 
for large e. We denote the two stable, stationary homogeneous states by {u±,v±), where 
the subscripts +(— ) refer to high (low) u and v values (see Fig. Ic), and ask what types 
of fronts connect these two states. 
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In the regime of large e, the field v varies on a time scale much shorter than that of 
■u, at least initially. Stretching the time coordinate by the factor e, we find 

Ut = v{u — — V -\- Uxx)i (3.6a) 

Vr = u — aiv — ao, (3.66) 

where u = <^ 1 and r = t/u. On that new scale, the leading order form of u is 
independent of r (set u = in (3.6a)) and equation (3.66) can be solved for v: 

v{x,t) = vo{x)e-^'''* + a^^[u{x,t) - ao][l - e"'"^*]. (3.7) 

Equation (3.7) implies that apart from a short transient of duration t ~ 0{e~^), the v field 
follows adiabatically the u field: v — {u — ao)/ai. Using this form in (3.6a) or (1.1a) we 
obtain the asymptotic system 

ut = (1 - a^^)u - - aoa^^ +Uxx, (3.8a) 

v={u — ao)/ai. {3.8b) 

Equation (3.8a) has the same variational structure as (3.1) has, and consequently, all the 
properties discussed in section 3.1 apply here as well. In particular, the property of fronts 
propagating in the direction dictated by the minimization of the free energy is retained. 
Note that the original system (1.1) is not variational, but reduces to a variational one in 
the large e regime. 

3.3 The regime < e < 1 

The regime of small e has been studied extensively (see the reviews [TyK88,Mer92]). 
The question we address here is how front solutions connecting the two states {u±,v±) 
differ from those we found for e » 1. 

Note first that for e <S 1 the v field remains approximately constant on the length scale 
over which u varies. Thus, in the narrow front region equation (3.1) applies with v = Vf, 
the value of v at the front position. The front speed is then determined by c = c{vf) 
where c{v) is given by (3.5). Imagine now that we prepare the system in the down state 
{u-,v-) and perturb it locally at the left edge of the system so as to induce a transition 
to the upper branch u = u+{v). This transition occurs at v = v_ and, as a result, a front 
propagating to the right, c = c(v_) > 0, will be induced as Fig. 4a demonstrates. If, on 
the other hand, the system is prepared in the up state {u+, v+) and perturbed on the right 
edge so as to induce a transition to the lower branch tt = tt_(v), a front propagating to 
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the left, c = c{v+) < 0, will be induced as shown in Fig. 46. Once the system is on a 
given branch it converges to the stationary homogeneous state lying on that branch. The 
two fronts therefore connect the same asymptotic states as x — > ±00, but propagate in 
opposite directions. This property of (1.1) reflects its nonvariational structure. Multiplicity 
of front solutions in reaction-diffusion systems with multiple time scales was first observed 
by Ortoleva and Ross [OrR75]. The stability of this type of solutions has been proved 
later by Rinzel and Terman [RiT82]. Note in this respect that the front solution of (3.8), 
calculated for e ^ 1, is independent of e and therefore continues to exist for e ^ 1 as well. 
In this regime, however, it is no longer stable. 

The picture revealed so far is as follows: For e » 1 there exists only one front solution 
having the asymptotic behavior {u,v) — > {u±,v±) as x — Too. For e <^ 1 three such front 
solutions coexist, two of which are stable. In the next section we study the intermediate e 
regime where we expect to find a bifurcation from a single to multiple front solutions. 

3.4 The front bifurcation 

We consider the symmetric case (ao = 0) for which front solutions of (3.8) are station- 
ary. In the following we will always assume the asymptotic behavior, {u, v) — > {u±, v±) as 
X =Foc), unless we specifically mention otherwise. The explicit form of the stationary 
front solution of (3.8) is (see (2.8)) 

Us{x) = —Uji- tanh.{r]x) , Vsi^) = ai^Us{x), (3-9) 

where = —u_ = — a^^), and ij = u+/\/2. We assume that at some e = €c 
propagating front solutions, u = Up{x), v = Vp{x), where x = — ct, bifurcate from 
(3.9). In the vicinity of the bifurcation point we can expand a propagating front solution 
in powers of its speed c: 

Upix) = Usix) + cuiix) + c^U2{x) + - 

Vpix) = Vsix) + cviix) + c^V2{x) + 

as well as e: 

e = ec + cei + c^e2 + .... (3-11) 
We insert these expansions into the propagating front equations 



'; + cu'p + Up-Vp-ul = Q ^^^^^ 
cVp + e{up - aiVp) — 0, 



collect all terms that contribute to a given order, and solve the resulting equations for 
the corrections ttj's and Vj's successively as described below. Note that the corrections to 
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all orders should decay to zero as x ^ ioo because the leading order form, that is the 
stationary solution, already satisfies the correct asymptotic limits. 

At order c we obtain: 

vi = a^^ui + u'g, (3.13) 

a^^ - Zul. (3.14) 

We consider a space of functions that decay to zero as x — ioo (e-g- Schwartz space). 
On this space M. is self adjoint and has the null vector u'^. 

Mu'^ = M^u'^ = 0. (3.15) 

Solvability of (3.13) then yields ec = aj"^. Using that result in (3.13a) we find A4ui = 
and thus ui = bUg, where b is arbitrary constant. Note, however, that u'g is a translation 
mode and that assigning a particular value to b defines the origin on the x axis. For 
simplicity we choose 6 = and thus obtain 

ui — 0, vi — u'g. (3.16) 

At order we obtain: 

= —eiofug + aiu'g, V2 = a^^'U2 + aiu". (3-17) 

Since Ug is odd, solvability of (3.17) requires ei = 0. Using that result in (3.17) we find 

U2 = yxK, V2 = ^xu's + Oi<- (3.18) 

At order we obtain: 

Mu3 = aju'l' - e2a?<, (3.19) 

which yields the solvability condition 62 = ~§(1 ~ '^i^)- 
Summing up these results we obtain 

Upix) =Us + \c^aixu'g + ... 

_ 1 (3.20) 

Vp{x) = ^1 + cu'g + (2^^'* + '^^'"v + 



where Ai is the linear operator 

M = Sld^ + 1 - 
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and 



e = ec + + 



-2 



e2 = -2(l-aJ-')/5. 



(3.21) 



We verified these results by comparing the pitchfork bifurcation imphed by (3.21) 
with a numerically computed one. The comparison is shown in Fig. 5a and indicates good 
agreement up to speed values c ~ 0{\^~^). Note that for the bistable case considered here 
ai > 1 and consequently ti is negative. This ensures that the bifurcation is supercritical. 

It is instructive to look at the propagating solutions (3.20) up to corrections of order 
c (including) . These can be written as 



Comparing with (3.9) we see that, to that order, the solutions for u and v remain the same 
except that v is translated with respect to u by amount proportional to the speed c. The 
direction of translation determines the direction of propagation as Figs. 6 demonstrate. 
This symmetry breaking is also reflected by the phase portraits of the front solutions in 
the u — V plane, as Fig. 7 shows. The stationary solution is described by a straight line, 
V — a^^u = 0, that does not break the odd symmetry of (1.1) (with ao = 0), while the 
propagating solutions break the odd symmetry and deviate from that line. This deviation 
is related to the front speed through 



where = 1/ J^^{u'^fdx. 

So far we have considered the symmetric model. In general the parameter ao will not 
be zero, unless there exists an inherent symmetry in the system that enforces that condi- 
tion. When ao 7^ the pitchfork front bifurcation unfolds into a saddle-node bifurcation 
as the bifurcation diagram shown in Fig. 56 suggests. The saddle-node bifurcation occurs 
at a critical value, ec(ao), smaller than Sc = a^^. 

3.5 Ising vs. Bloch fronts 

An interesting distinction between the stationary and the propagating front solutions 
in the symmetric model can be made once we define an angle or phase in the u — v plane, 
9? = arctan(w/w). Across the stationary front the phase is constant except at the core 
where it suffers a jump by tt and the norm {v? + f^)^/^ vanishes. Across a propagating 
front, on the other hand, the phase rotates smoothly by an angle tt and the norm never 
vanishes. An analogous front bifurcation has been found recently by CouUet et al. [CLH90] 
in the forced complex Ginzburg-Landau (CGL) equation 



«p(x) ~ '"s(x), 



^p(x) ~ ^«s(x + cai). 



(3.22) 




At = {p + id)A + (1 + ia)A^^ - (1 + i(3)\A\''A + 7^, 



(3.23) 
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where A{x,t) is a complex field. This equation describes spatiotemporal modulations of 
an oscillating medium, periodically forced at twice the oscillation frequency with strength 
7. The effect of the forcing is to create two stable, phased locked states whose phases, or 
arg(^) values, differ by tt. The phase arg(^) plays the same role as (p does in the FHN 
model; it remains constant or smoothly rotates across a front separating different phase 
locked states for sufficiently large or small 7 values, respectively. CouUet et al. [CLH90] 
proposed these types of fronts solutions to be the nonequilibrium analogs of Ising and 
Bloch walls in ferromagnets with weak anisotropy (where arg(^) corresponds to the angle 
the magnetization vector makes with an easy magnetization direction) [LaN79]. In the 
following we will occasionally follow this terminology and refer to the two types of fronts 
in the FHN model as Ising and Bloch fronts. 

3.6 The bifurcation line in the e — 6 plane 

The extension of the bifurcation analysis of section 3.4 to nonzero d values is rather 
tedious and we evaluate here the bifurcation line only. Even this calculation is not simple 
as the exact stationary front solution is not known. We therefore consider the regime of 
small 6 values where the stationary front solution can be expanded in powers of S, and 
evaluate a linear approximation for the bifurcation line valid for 5 -C 1. At the other 
extreme, d/e 3> 1, we use a singular perturbation analysis, similar to that of Ikeda et al. 
[IMN89]. The complete bifurcation line is obtained numerically by integrating (1.1) with 
Neumann boundary conditions. 

Consider the symmetric model, that is, (1.1) with ao = 0. A stationary front solution 
u = Us{x; S) V = Vs{x;S) satisfies 



Note that in these notations u = Us{x; 0), v = Vs{x; 0) is the stationary solution (3.9). For 
5 <^ 1 we can solve (3.24) perturbatively. Since we shall need the form of the stationary 
solution on the bifurcation line we set e = ec{5) in (3.246). Expanding the stationary 
solution, u = Us{x; S) v = Vs{x; S), and ec{S) in powers of S, 



+ Us - Vs - ul = 0, 



(3.24a) 



Sv" + e{us — aiVs) = 0. 



(3.246) 



Us{x; 6) = Us{x; 0) + 6ui{x) + .. 
Vs(x; S) = Vs(x; 0) + Svi(x) + .. 



(3.25) 



ec{S) = oi ^ + Ssci + 



(3.26) 
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and using these expansions in (3.24) we find at order 6 

Mui{x) = u'^{x;0), (3.27) 

vi{x) — ai^ui{x) + Ug{x;0), (3.28) 

where M. is given by (3.14) with Ug = Us{x; 0). The solution of (3.27) is ui{x) = ^xu'g{x; 0). 
This, together with (3.28), yields the leading order approximation 

1 

Ug{x; 5) — Us{x; 0) + -5xu'g{x; 0) + (3.29a) 

Vs{x\ 8) = a^^Usix; 0) + ^5a^^xu'^{x; 0) + 5u"{x; 0) + (3.296) 

where tts(a;;0) is given by (3.9). 

As before we expand e and the propagation front solutions in power series in c (see 
(3.10)) 

«p(x; = ^slx; + cwi(x; <^) + 
^p(x; = ^s{x\ ^) + cvi(x; S) + 

e = ee(5) +cei(5) + ..., (3.31) 
and insert these expansions into (1.1). At order c we obtain 

^2 I 1 Q„,2 



(3.30) 



a^ + l-3«^(x;5) -1 \ u^{x\^) 



ed 5 d^- ecai/5 J \ vi (x; S) 



(3.32) 



Solvability of (3.32) requires the right hand side of this equation to be orthogonal to the 
vector {u'g{x'i ~j~'"six'i ^))^- This yields the following equation for the bifurcation line: 

/oo t-oo 
v's{x;5fdx u',{x;Sfdx. (3.33) 

-oo J — oo 

Using the leading order approximation (3.29) in (3.33) we find 

e^{6) = {l-bS)/al or ^^(e) = (1 - a?e)/6 {S<^1), (3.34) 

where b = 4(ai — l)/5. 

In Appendix A we evaluate the bifurcation line for 5/e ^ 1 using singular perturbation 
theory. This approach differs from that presented above in that it uses the smallness of 
e/5 rather than that of the speed. As noted by Ikeda et al. [IMN89] it yields approximate 
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front solutions and bifurcation diagrams valid away from the bifurcation point as well. The 
price though is the limited validity range along the bifurcation line. Using this approach 
we obtain 

ec(S) = d/5, d = 9/(2ai + 1)^ (3.35) 

where we took ao = (symmetric model). 

Fig. 8 shows the exact bifurcation line obtained by direct numerical integration of 
(1.1). The linear approximation for 5^1 and the singular perturbation result (3.35) 
(valid for 6/e ^ 1) are shown in solid and dashed lines, respectively. We conclude this 
section by showing in Figs. 9 typical phase portraits of the front solutions (more precisely, 
projections thereof on the (tt, v) plane) for 6 values of order unity or larger. The outer 
portraits pertain to the two counter propagating front solutions that exist for e < ec(5), 
whereas the middle one represents the stationary front solution at e > €c{6). Note that 
the field v barely changes across the front. This is a consequence of the smallness of e/5 
(compare with the corresponding phase portraits shown in Fig. 7 for 5 = 0). 

4. Implications on Pattern Formation 

4.1 The Emergence of Persistent Patterns 

A significant implication of the front bifurcation for 5 = or sufficiently small is that 
it tells us where in parameter space we should expect an initial pattern of domains to 
decay toward a uniform state and where to develop into a stable traveling wave. Consider 
first the regime of high e values, e ^ ec(ao) (see last paragraph of section 3.4). In this 
parameter regime the nonvariational system (1.1) reduces to the variational system (3.8). 
As a result, patterns evolve so as to minimize the appropriate free energy functional. Like 
in equilibrium phase transitions, up state domains either shrink or expand (depending on 
the sign of ao) to form the uniform state of lowest free energy. The relaxation toward 
a uniform state occurs even for the symmetric case (ao = 0) where an isolated front is 
stationary. As shown by Kawasaki and Ohta [Ka082] and more recently by Carr and Pego 
[CaP89] , nearby fronts that bound a domain attract and annihilate one another on a time 
scale t ~ exp{X/w) where A is the domain size and w is the front width. Thus, domains 
shrink and disappear, but the relaxation can be extremely slow and unnoticible in practice. 

We have found numerically that this qualitative behavior remains unchanged for all 
e > ec(ao). Fig. 10a shows the relaxation of an initial pattern toward a uniform state 
for e/ec(ao) = 1.12. Fig. 10b shows the time evolution of the same initial pattern for 
e/ec(ao) = .72, that is beyond the bifurcation from Ising to Bloch fronts. In contrast 
to Fig. 10a, here, a stable traveling pattern develops. We explored the transition from 
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transient to persistent patterns more carefully by integrating equations (1.1) for various e 
and ao values keeping the other parameters constant {5 = 0, ai = 2). The resulting phase 
diagram is shown in Fig. 11. The solid (right) curve designates the front bifurcation, 
e = ec(ao). To the right of this curve only one front solution {u, v) — > {u±, v±) as a; — > ^oo 
exists, and initial patterns decay toward a uniform state. To the left of this curve, three 
front solutions connecting the same asymptotic states coexist; two stable (Bloch fronts) 
and one unstable (Ising front). Stable traveling waves, however, appear only beyond a 
second threshold, ep(ao), given by the dashed (left) curve. Using a numerical continuation 
method we found that for the symmetric case (ao = 0) traveling wave solutions do exist for 
ep < e < ec, but they are unstable. We did not find such solutions for the nonsymmetric 
model (i.e. for ep(ao) < e < ec(ao) with ao ^ 0). 

This dramatic change in the qualitative behavior of patterns can be attributed to 
two related factors; the appearance of front multiplicity and the appearance of a second 
independent field. The former makes it possible for domains to travel rather than shrink 
or expand, whereas the latter affects front interactions so as to bind a trailing front to a 
leading one. In the rest of this section we elaborate on these two factors in some length. 

The multiplicity of front solutions for e < ec(ao) and the symmetry x — > —x of (1.1) 
imply that along with a front that transforms the lower state (■u_,v_) to the upper state 
{u+,v+), there exists another front (hereafter "back") propagating in the same direction 
that transforms the upper state back to the lower one (note that the front and the back 
represent two difi^erent Bloch fronts). A combination of the two, using appropriate initial 
conditions, can yield a traveling up-state domain. To be concrete, let us assume that the 
(isolated) trailing front, or back, is faster than the leading front. This situation is attenable 
with ao < 0. The fate of the traveling domain as the back approaches the (leading) front 
depends on the interaction between the two. Unlike the case of Ising fronts where the field 
v{x, t) is eliminated through the relation v = (u — ao)/ai and the interaction is attractive, 
for e < ec(ao), v{x,t) is independent of u{x,t) and allows for a repulsive interaction as 
well. 

To see this, let us assume that e <^ Sp where < ec ~ 0{a^^). Then v varies on time 
and length scales much longer than those of u. Just behind the front the value of v is still 
close to the down state value It approches the up state value, v+, only far behind the 
front. Since v is approximately constant across the narrow back region, the back speed is 
determined by the local v value, Vb, according to (3.5). When the back is still far behind 
the front Vf, v+, the energy difference V(tt_(v+)) — £{u+{v+)) is large (see Fig. 3c) and 
the speed high. As it approaches the front, however, Vh decreases, the energy difference 
becomes smaller and the speed lower. Thus, the back is slowed down (or repelled) by the 
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field V induced by the front. For e sufficiently small the back approaches the front speed 
while it is still far behind the front. At such distances the u field follows adiabatically 
the slow field v (that is, u = u^{v)) and a stable traveling domain of fixed size is formed 
(see below). For e values closer to ec the adiabatic elimination of u may not be valid and 
more complicated interactions are possible. We believe that in this parameter range both 
repulsive and attractive interactions are significant and suggest that the balance between 
the two gives the threshold ep represented in Fig. 11 by the inner curve. We note though 
that we have not studied yet this regime in detail and that the above suggestion is not 
founded yet. 

To study the formation of a stable traveling domain of fixed size we assume e <^ e^, 
and use a singular perturbation approach. We first evaluate the slow field v that builds 
up behind a front propagating at constant speed c = c(v_). In this region u = u+{v) and 
the equation we need to solve is 

cv' + e(u+{v) — aiv — ao) = 0, (4-1) 

where v = v{x) and x = a; — ct. The boundary conditions are 

v{x^f)^v-, v{x)^v+ as x^-oo, (4.2) 

where %/ is the front position in a frame moving at speed c. To simplify (4.1) we take ai 
to be sufficiently large so that \v- \ ~ v+ <S 1. The branches u± can then be linearized: 

u±{v) = ±l-v/2. (4.3) 

Using (4.3) in (4.1) and solving for v we obtain 

v{x) = {v. - v+)e^<^-^"f^ +V+, x< X/, (4.4a) 

where 

±l-oo ai + 1/2 

v± = — —7-, K = . (4.46) 

ai + 1/2 c 

A traveling domain of fixed size can be obtained by demanding the front and back 

speeds to be equal 

-c(w°) = c(w_) = c. (4.5) 

This relation together with (3.5) can be used to evaluate = ^(Xb)) the level of v at 
the back position for a domain of fixed size. The fixed size of an up-state domain can be 
expressed in terms of v° using (4.4a): 

A = X?-xg = -ln(^^^^). (4.6) 
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Note that (4.6) requires v° < v+ or an isolated back to propagate faster than an isolated 
front. For v° > v+ (back slower than front) the v field behind the front can never reach the 
value v° needed for a back to propagate at the front speed. In this case, an up-state domain 
expands indefinitely, and a down-state domain of fixed size (referred to as an E-pulse in 
[RiT82]) becomes feasable. We illustrate these two cases in Figs. 12. 

The relaxation toward an up-state domain of fixed size can be derived by writing a 
back solution in the form 

Ub{x,t) =u'l{x-Xb) + eul{x-Xb,et), (4.7a) 
Vb = v{xb) = v^ + evl{et), (4.76) 



where w^(x) solves 



uf + cul' + ul- vl - uf = 0, (4.8a) 
«g(x) = (^-!l| ""^'^ (4.86) 



and 



Xh{et) = xl + Xb{et) 



is the actual back position. The slow time dependence of Vh is introduced to account for 
the declining v field at the back position as the latter approaches the front. Using (4.7) in 
(1.1a) we find 

dlul + cd^, + (1 - Zuf)ul = vl - e-'xb^, (4.9) 

where the dot over Xb denotes differentiation with respect to the fast time t. Solvability of 
(4.9) requires the right hand side of to be orthogonal to exp{cx)du1/dx or 



Xb = Pevl (4.10) 

where 

is a positive constant. 

Equation (4.10) relates the back speed to the level of v at the back position. The 
latter, in turn, can be related to former through (4.4a): 

evl = v{Xb) - v'b = {v+ - ^6)(1 - e'^^"). (4.11) 
Combining (4.10) and (4.11) we obtain an equation of motion for the back: 

ib = Piv+-v^b){'^-e'^^')- (4.12) 
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The linearization of (4.12) about X6 = (i.e. the fixed size domain) gives the equation 

Xb = -eK(3{v+ - vl)xh- 

Recall that e, and (3 are all positive and < vj^ for a relaxation toward an up-state 
domain. We therfore conclude that the up-state domain whose size is given by (4.6) is 
stable to translational perturbations. Note that in the small e regime considered here the 
size of a domain is much larger than the front (back) width. Thus, the front and back that 
bound a domain are approximately isolated (see (4.7a). The isolated front solutions, in 
turn, are marginally stable to translations but stable to other types of perturbations (such 
as amplitude modulations). We therefore expect the translational perturbations to be the 
most "dangerous" ones. 

A similar approach can be used to study the relaxation of arrays of traveling up- 
state domains toward periodic traveling waves. The difference with respect to a single 
domain is that now the leading front of any domain propagates into the slowly decaying 
V field behind the back ahead of it. As a result, the front speed is not constant; it rather 
depends on the distance between the front and the leading back. A different kinematical 
approach, leading to equations of motion expressed in implicit form, has been proposed 
earlier by Keener [KeeSO] in the context of excitable media. We also note that the singular 
perturbation theory sketched above cannot capture an oscillatory decay toward either of 
the stationary homogeneous states. Such a decay, can lead to spatial chaos as discussed in 
refs. [EMS88,Mer92]. 

So far we have considered the case 5 = 0. Similar conclusions hold for small enough 
5 values. We conclude this section with a numerical study of the forced CGL equation 
(3.23) for which 5 = 1. We have added to the right hand side of this equation a constant 
term, 70, to break the symmetry between the two phase locked states (analogous to ao). 
Such a term appears when the periodic forcing contains a component at the medium's 
oscillation frequency. We have not explored in detail the (7, 70) parameter plane, but for 
fixed 7o we did observe transient patterns for 7 sufficiently large, stable traveling waves for 
7 sufficiently small, and transient patterns involving traveling domains in a narrow range 
of intermediate 7 values. This suggests that the qualitative form of the phase diagram 
drawn in Fig. 11 is quite general and might apply to any system (with 5 relatively small) 
undergoing an Ising-Bloch type front bifurcation. 

4.2 Stationary, Oscillating and Traveling Domains in an Infinite System 

When the diffusion of v is fast enough stable stationary domains, rather than transient 
ones, may develop [EHT84] in the regime of Ising fronts. The formation of this type of 
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domain for 5/e sufficiently large is illustrated in Fig. 13. As the two fronts approach one 
another the v field, inside the domain bounded by the fronts, developes a narrow lumped 
structure. The consequent diffusive damping lowers the level of v at the front regions and 
thus slow them down. The size. A, of a single stationary domain is given by 

A = w{6/ey/^ w=- ln(-ao)/(ai + ^f'^ (4.13) 

assuming fx = e/S <^ 1. We derive this result in Appendix B using a singular perturbation 
approach. Note that the stationary domain size decreases as iJ,~^ = 5/e decreases. Obvi- 
ously, for such a domain to exist its size should be at least of the order of the front width. 
This consideration explains the existence of a critical value /ijf^, below which stationary 
domains do not exist. We note that the transition to stationary patterns beyond //^t^ does 
not occur by a Turing bifurcation. We remind the reader that we are considering here case 
{in) pertaining to bistable systems, and that none of the homogeneous states of the system 
lose stability at the transition point. The transition line in the e — 5 plane, 5 = iJ,J^e, does 
resemble, however, the Turing bifurcation line (2.1). 

To study how patterns are affected by the front bifurcation, we begin with a single 
stationary domain in an infinite medium and follow that structure as we approach the 
front bifurcation line. We postpone the discussion of periodic domain patterns to the next 
section. The results to be described below have been obtained numerically. For more 
details about the numerical procedures the reader is referred to section 5. Our starting 
point in the e — 6 plane lies above the transition to stationary domains, 5/e > fij^, and 
far to the right of the front bifurcation line, namely, at e values significantly larger than 
ec{5) (see (3.34) and (3.35)). As e is decreased (keeping 5 constant) past a critical value, 
ebri^), the stationary domain loses stability in a Hopf bifurcation, and a breathing like 
oscillatory motion sets in as shown in Fig. 14a. This type of oscillations have first been 
observed by Koga and Kuramoto [K0K8O] and more recently have been studied analytically 
by Nishiura and Mimura [NiM89] (see also [KeL90], [OIT90]). In the latter study a finite 
medium has been considered and the oscillations were found to persist well below ebr 
without any indication of a secondary bifurcation. Our findings for a single domain in 
an infinite medium are different; as e is further decreased the oscillation amplitude grows, 
leading eventually to the collapse of the domain and to a uniform down state as Fig. 146 
illustrates (note that the homogeneous states {u±,v±) remain linearly stable for any e and 
S values). Similar collapse events have recently been found in a model for a semiconductor 
etalon [RRH93]. The collapse can be avoided by increasing the size of the domain. This 
can be achieved by making the system more symmetric, that is, by decreasing |ao|. For 
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|ao| sufficiently small oscillating domains have been observed all the way down to the front 
bifurcation line, e = e^S), but not below it! 

The effect of the front bifurcation on the behavior of a single domain structure is 
demonstrated in Figs. 15a and 156 corresponding to e > ec{5) (Ising fronts) and e < ec{5) 
(Bloch fronts), respectively. In both figures the initial conditions consist of two fronts 
bounding a wide up-state domain and propagating toward one another. In 15a the two 
fronts set in a stable oscillatory motion, whereas in 15b they rebound from one another 
and leave the system in a uniform up state. We attribute this change of behavior to the 
appearance of multiple fronts. Above the bifurcation line (e > ec{S)) there exists only one 
type of stable front solution, approaching {u±,v±) as x — > ^oo and propagating to the 
left (assuming ao negative), or, the symmetric one, approaching (tt=p,v=p) as x — > ^oo and 
propagating to the right. Below the bifurcation (e < ec{5)), another pair of stable front 
solutions appear (in addition to the aforementioned one): a front approaching {u±,v±) 
^ X ^ -poo cind propagating to the right, and a front approaching (tt=p,v=p) as x 
and propagating to the left. In 15b the outward propagating front structures converge 
to the second pair of stable front solutions and thus proceed with the outward motion 
indefinitely. In 15a, on the other hand, there exist no front solutions that pertain to 
an outward motion. The only front solutions that exist are those propagating toward 
one another, and the change in the direction of propagation just reflects the convergence 
toward these solutions. 

Figs. 16 show another numerical experiment where the same parameter settings are 
used but the initial conditions consist now of two fronts propagating in the same direction. 
For e > ec{5) (Fig. 16a), the front on the right side of the domain changes its direction of 
propagation and an oscillatory motion sets in. For e < ec{5) (Fig. 166), a stable traveling 
domain, is obtained. 

We note that the stationary domain solution still exists on both sides of the bifurcation 
line but is unstable. Fig. 17 shows a typical bifurcation diagram for domain structures. It 
does not capture the oscillatory domain solutions or the collapsing solutions (see section 
5) but it shows the destabilization of the stationary domain solution at e = e^r, and how 
traveling domain solutions bifurcate from it. We also note that by making the system 
more symmetric (i.e. decreasing |ao|) the destabilization of the stationary domain solution 
occurs closer to the bifurcation to traveling domains. For sufficiently symmetric systems 
it might be possible for the stationary domain solution to remain stable all the way down 
to the traveling-domain branch, but we have not verified that. 

A bifurcation from stationary to traveling domain solutions has already been found in 
the context of excitable media [DoK89] . What is new in our findings is that this transition 
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is a consequence of the front bifurcation. Indeed, we have verified numerically that the 
appearance of traveling domain solutions at e = etr{S) (see Fig. 17) coincides, within the 
numerical accuracy, with the front bifurcation line, e = edS), at least above the onset of 
stationary domains (the line 5 = /J^Jt^e). Well below that line etr deviates from Sc and 
converges, at 5 = 0, to (see section 4.1). 

We summarize the discussion of single domain structures in infinite media by showing 
in Fig. 18 a phase diagram in the e — S plane. The diagram is divided into four main 
regions. Region I, e > etr{S) and e > jist^, where domain structures are transient and the 
asymptotic state is either the uniform up or uniform down state. Region II, e > €5^(5) and 
e < jist^j where stationary domains prevails. Region III, e > etr{S), e < fJ,st^j ^ind 
e < ebr{S), where domains oscillate. And region IV, e < etr{S), where domains travel 
invariably. The hatched area in region III designates, in a schematic manner, the regime 
of steadily oscillating domains. Outside this regime domains oscillate for a while and then 
collapse. Depending on the asymmetry of the system, this regime may shrink to a narrow 
strip below the line e = ebri^) (strong asymmetry), or span the whole range between the 
lines e = e^riS) and e = etr{S) (weak asymmetry). The range between the two lines, in 
turn, becomes smaller the more symmetric is the system. Note the thick solid-dashed line. 
This is the front bifurcation line (e = ec{5)) which coincides with the onset of traveling 
domains at e = etriS) apart from a small portion at small S values. 

4.3 Periodic Domain Patterns 

The behavior of periodic domain patterns is similar in many respects to that of single 
domain structures. One notable difference though is the persistent of stable oscillating 
domain solutions to the left of the front bifurcation line (e < ec{5)), where they coexist 
with stable traveling domain solutions as Figs. 19 demonstrate. The regime of coexistence 
seems to be limited; far enough from the front bifurcation line we found the traveling 
domain solutions to prevail. The existence of oscillatory dynamics in the regime of multiple 
fronts is a consequence of front interactions. As nearby domains expand and approach one 
another there is little room left for diffusion of v. As a result the level of v at any front 
position increases and the fronts change their directions of propagation. 

Another difference between periodic and single domain patterns is the possibility of 
having various modes of oscillation. Ohta et al. [OIT90] have identified three such modes; 
an in-phase mode where all domains expand and shrink simultaneously, an out-of-phase 
mode where nearby domains oscillate with a phase shift of tt, and an "acoustic" mode, 
where nearby domains travel in opposite directions back and forth (note that in the latter 
case the spaces between the domains, or the down-state domains, undergo out-of-phase 
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oscillations). Of the three, the in-phase mode is the first to be excited. As pointed out 
by Ohta et al. [OIT90], this can be understood by considering the stabihzing effect the 
diffusion of v has. The in-phase mode involves uniform deviations of v from the stationary 
domain profile and stabilization by diffusion does not take place. The out-of-phase mode, 
on the other hand, pertains to alternating deviations every second front site. In this case 
diffusion of v tends to smear out these deviations and thus to stabilize the stationary 
pattern against this mode. Similar considerations hold for the acoustic mode. In the 
numerical experiments we have conducted the in-phase mode always prevailed, even at 
relatively small 5 values, where the diffusive damping of that mode is weaker. 

We expect Hopf-Turing systems, driven far beyond the Turing instability, to exhibit 
the same phenomenology that short-wavelength patterns in bistable systems do. Indeed, 
by decreasing e at constant S (far beyond the Turing instability) we found a transition 
from stationary to oscillating patterns, followed by a transition to traveling patterns. 

4.4 Reflection at Boundaries 

The front bifurcation may drastically affect the behavior of fronts near an impermeable 
boundary. The effect of such a boundary is similar to that of another approaching front 
and the scenario of behaviors is similar to that found for domain structures. Figs. 20 
show the boundary effects as we cross all four regions in the e — S plane by decreasing e 
at constant 6 (the simulations were carried out with Neumann boundary conditions). In 
region I (transient domains) the front is absorbed at the boundary, leaving the system 
in a uniform state. In region II (stationary domains), the front comes to a stop at a 
characteristic distance from the boundary, whereas in region III (oscillating domains) it 
oscillates near the boundary. In all these regions the (Ising) front is either absorbed at 
the boundary or is bound to it. By contrast, in region IV the front is reflected. This is a 
consequence of the coexistence of two counter propagating Bloch fronts beyond the front 
bifurcation line. At front speeds high enough relative to the diffusion of v, both Ising 
and Bloch fronts are absorbed at the boundaries. Ising fronts, however, never reflect at 
boundaries. 

5. Numerical Procedures 

Numerical simulations of (1.1) were performed using the method of lines with the 
spatial derivatives approximated by 2nd or 4th order finite differences on a uniform mesh. 
To solve the resulting system of ordinary differential equations, we used a stiff ODE solver 
[NKM89]. This solver implements the implicit methods of Gear [GeaTl] with adaptive 
control of both time step and method order. 

For the computations pertaining to single domain structures in inifinite systems the 
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simulations were done with Neumann (no flux) boundary conditions and large computa- 
tional system sizes. We checked the results with larger system sizes to ensure that the 
boundaries did not contribute any measurable effect on the solution. When the simula- 
tions pertain to periodic patterns, periodic boundary conditions were applied. Occasionally 
we used periodic boundary conditions for single traveling domains rather than Neumann 
boundary conditions (see for example Fig. 166). In those cases we chose the period to 
be large in comparison with the domain size and checked that the two types of boundary 
conditions gave the same solution within the numerical accuracy. 

For the continuation of solutions we used the bifurcation package AUTO [Doe86]. 
Starting with a stable traveling or stationary solution computed with our PDE solver we 
continued a periodic solution in one or two parameters to study the effects of parameter 
variations on the existence and stability of solutions. When AUTO failed to give accurate 
stability information, due to the approach of the periodic orbit towards two heteroclinic 
cycles, we used our PDE solver to integrate the computed solution forward in time to check 
its stability. 

To find the oscillating domains region we linearized (1.1) about the numerically com- 
puted stationary solution and solved the resulting eigenvalue problem. The spatial part of 
the linearized operator was discretized using 2nd order finite differences on a nonuniform 
mesh with grid points concentrated near the steep gradient front structures of the domain. 
We solved the eigenvalue problem using the IMSL routine DEVCRG for real general matri- 
ces. Decreasing e below the line e = ebr{S) a pair of complex conjugate eigenvalues crosses 
the imaginary axis indicating that the stationary stable solution goes unstable through a 
Hopf bifurcation to an oscillating domain. 

6. Conclusion 

We have studied a front bifurcation, likely to occur in a variety of nonvariational 
reaction diffusion systems, and the possible consequences of this bifurcation on pattern 
formation and pattern dynamics. The bifurcation is accompanied by the appearance of 
a second independent field which breaks the odd symmetry of front solutions, and leads 
to a richer dynamical behavior. The main effect of the bifurcation, regarding patterns, is 
to allow the formation of stable traveling domain patterns. For small 5 values (including 
5 = 0), it provides a mechanism for pattern formation as no persistent domain patterns 
can be formed below the bifurcation. For larger 5 values it is responsible for transitions 
from transient, oscillating, or stationary patterns to traveling ones. Other consequences 
of the bifurcation are a coexistence regime of traveling and oscillating patterns, and the 
possible reflection of fronts at boundaries. 



- 24- 



These results were obtained for bistable systems but are also relevant to excitable 
and Hopf- Turing systems as the latter can be obtained by unfolding the former. In this 
sense, the variety of chemical-wave phenomena in excitable and oscillatory systems can be 
regarded as a consequence of an Ising-Bloch front bifurcation occurring nearby in parameter 
space. A more specific implication for Hopf- Turing systems is the prediction of a transition 
from stationary to oscillating patterns, and from oscillating to traveling patterns far beyond 
the Turing instability. These transitions may shed light on the complex spatio-temporal 
behavior that has recently been observed in the CIMA reaction by Ouyang and Swinney 
[OuS91]. 

Beyond these specific results and implications, the present analysis demonstrates 
the usefulness of focusing attention on localized structures that appear in extended pat- 
terns. This view lies behind the various kinematical approaches to pattern dynamics 
[KiM89,NeM92,Mer92,CrH93], where the effects of small perturbations and weak interac- 
tions on the dynamics of localized structures are studied. The point we wish to emphasize 
here is that it is also important to look for structural changes that localized structures 
might undergo as these can drastically affect pattern dynamics. 

We confined ourselves in this paper to patterns in one space dimension, but we expect 
the front bifurcation to have important effects in two dimensions as well. Perhaps the first 
question that needs to be addressed concerns the dynamics of a curved Bloch front near 
the front bifurcation point. The dynamics of a such a front far beyond the bifurcation is 
well known [TyK88,Mer92]; the normal speed is affected by the curvature k according 
to the simple linear law, = c — Dk, where Z) is a positive constant having the dimension 
of a diffusion coefficient. The consequent effect of curvature is to smooth out any wiggles 
along the front and thus to stabilize its dynamics. Close to the bifurcation point we expect 
this relation to be highly nonlinear. This might change the role of curvature in a way that 
still needs to be explored. 

An important consequence of the Ising-Bloch front bifurcation in two space dimensions 
is the possible formation of spiral waves. As pointed out by CouUet, a spiral wave can 
be created by joining end to end two counter propagating Bloch fronts. This view goes 
back to Fife [FifSS] who discussed the twisting action induced by joining a front and a 
back end to end in an excitable medium. All studies of spiral waves so far have taken e 
to be a small parameter, and therefore are only valid far from the front bifurcation. The 
expected nonlinear relation between curvature and normal speed close to the bifurcation 
point makes the study of spiral waves in that regime potentially interesting. 
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Appendix A: The front bifurcation for e/(5 <S 1 

We consider here the regime t/5 and use singular perturbation theory to derive 
the front bifurcation hne. Rescahng space and time according to 

z — \f]lx^ T = et, // = e/5-Cl, 

and assuming travehng solutions, u = u{() and v = v{(), where = z — ct, we obtain 
from (1.1) 

uucc + cSuuc + u — u^— v = 

{A2) 

V(^(^ + CV(^ + u — aiv — oo = 0. 

We consider front solutions of {A2) that satisfy (ti, v) {u±, v±) as C — * =Foo, and set the 
origin, C = 0, by demanding u{0) = 0. Outer solutions of {A2) are obtained by setting 
H = and solving 

V(^(^ + CV(^ + u±{v) — aiv — ao = 0, (^3) 

for C < (tt = u+{v)) and ( > {u = U-{v)). To simplify (^3) we approximate the 
branches u±{v) by the linear forms (4.3). This approximation can be controlled by varying 
ai. We then obtain the following boundary value problems: 

^ + ' s ^^^^ 

v{ — Oo) = V+ v{0) = Vf, 

and 

^ (^5) 
f (oo) = V- f (0) = f/, 

where v± are given by (4.46), A;^ = ai + | and is the level of v at the front position. 
The solutions are 

v{C) = {vf-v+)e'''^ + v+, C<0, (^6) 

^(C) = (^/-^^-)e"^'^ + ^-, C>0, {A7) 

where 

ai,2 = -c/2±(cV4 + fc')'/^ (^8) 
Matching the derivatives of the two outer solutions at ^ = yields 

c ao 
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A second relation between Vf and c is obtained by solving the inner problem. To this end 
we stretch the traveling-frame coordinate according to x = C/v^ obtain from (^2) 

Uyy + VCU-y + u — u — V = 0, 

(AlO) 

where r]^ — eS. Setting formally /j, — in {AlO) we obtain a nonlinear eigenvalue problem 
for c, 

u^^ + rjcu^ + u — u — Vf = 0, 

(All) 

•u(=Fcxd) = u±{vf) = ±1 — Vf/2, 
which yields the second relation between Vf and c: 

/2 

Vf = -^Vc. {A12) 

Comparing {A9) and (^12) we find 

3^ 2A;2(cV4 + fc2)i/2 + ^^^^^ 

A plot of the solutions c = c{r]) of (A13) in the (c, rj) plane yields a bifurcation diagram with 
rj being the bifurcation parameter. Consider now the symmetric case, ao = 0. Assuming 
propagating solutions (c ^ 0) and taking the limit c ^ we find the bifurcation point 



Vc = 



(2ai + 1)3/2' 
or, recalling that rj^ = eS, the bifurcation line 
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Appendix B: Stationary-Domain Solutions for e/5<^l 
Stationary-domain solutions of (1) satisfy 

UUzz + u — u^— v = ^ 

{B2) 

Vzz + u- aiv - ao = 0, 

where z = ^JJix and = ejb is sufficiently small. In the following we will assume that 
<^ 1 and that ao < 0. The latter assumption means that an Ising front connecting an 
up state on the left with a down state on the right propagates to the left and consequently 
that a single domain is an up-state domain (see Fig. 13). Outer solutions are obtained 
by setting = in (-B2). We look for such solutions in three regions, z < zi^ z > z^^ 
and zi < z < z^^ where zi and z^ are the positions of the fronts that bound the stationary 
domain from the left and the right, respectively. The first two regions pertain to down-state 
domains in which we need to solve 

Vzz - k'^v + k'^v^ = 0, {B2) 

subject to the boundary conditions 

v(—oo) = V- v(zi) = vi for z < zi, 

(S3) 

v{oo) = f_ z{Zr) = Vr for Z > Zr, 

where vi and Vr are yet to be determined. The third region, zi < z < Zr, corresponds to 
an up-state domain in which we solve 

- k\ + k'^v+ = 0, (S4) 



subject to the boundary conditions 

v{zi) = Vi v{Zr) = Vr- (-B5) 

In {B2) and (-B4), /c^ = ai + 1/2 and we have used the linearized forms of the branches 
u±{v) (see (4.3)). 

The solutions of (S2) to (S5) are: 

v{z) = {vi - v-)e^^^-^'^ + v_ z< zi, {B6) 



v{z) = (vr - v-)e ^''^ ^'■^ -h V- 
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and 

v{z) = . ,\ (vr -v+) smb. k{z - zi) - {vi -v+) smb. k{z - Zr) + v+ ^ ^ 
smh(/cA) I V / V -r/ 'J ^ (^^8^ 

Zi < Z < Zri 

where K = Zr — zi. Matching the first derivatives of these solutions at z = zi and z = Zr 
we find that vi = and that A (the domain width) solves the equation 

coth /cA — csch/cA = {vr — V-)/{v+ — Vr). 

Equation (-B9) still contains one unknown parameter, Vr- This parameter is obtained 
from the analysis of the inner regions centered at ^ = and z = Zr- Stretching the 
coordinate system according to a; = zj -^f]! and setting = we find 

'^xx -\- u — V? — Vr = ^1 (-^10) 
implying Vr = 0. Using this result in (S9) and solving for A we get 

A=-^ln(-ao), (Sll) 

or (4.13) where A = h./ yfjl. 
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Figure Captions 

Figure 1: Three basic cases for the nuUchnes of the model system (1.1). Case (i): The 
nuUchnes intersect at a single point on an outer branch of the cubic nuUcline. Case (ii): 
The nuUclines intersect at a single point on the middle branch of the cubic nuUcline. Case 
(iii): The nuUclines intersect at three points each on a different branch of the cubic nuU- 
cline. 

Figure 2: Patterns near and far from onset show the development of two spatial scales, 
(a) Stationary periodic pattern near onset of the Turing instability, (b) Stationary peri- 
odic pattern far from onset of the Turing instability, (c) Traveling pattern near the Hopf 
bifurcation, (d) Traveling pattern far from the Hopf bifurcation. In all the figures the solid 
line is the u field and the dashed line is the v field. 

Figure 3: The free energy density £ (see (3.3)) for different values of v. 

Figure 4: Multiplicity of stable front solutions when e <^ 1. (a) Preparing the system 

in the down state and perturbing it locally at the left edge of the domain induces a front 

propagating to the right, (b) Preparing the system in the up state and perturbing it locally 

at the right edge of the domain induces a front propagating to the left. 

Figure 5: Bifurcation diagrams of front solutions. The dots are data points representing 

the speed of the different types of stable front solutions that exist for each value of e. 

(a) The symmetric case (ao = 0). The solid line is the theoretical bifurcation diagram 

computed from (3.21). (b) The nonsymmetric case (ao = .1) shows the unfolding of the 

pitchfork to a saddle node bifurcation. 

Figure 6: The propogation direction of Bloch fronts is determined by the translation of 
the V field relative to the u field; v always lags behind u unless the front is stationary, (a) 
Left traveling front, (b) Stationary front, (c) Right traveling front. 
Figure 7: Phase portraits of front solutions connecting the {u+,v+) state at x = — 
to the (■u_,v_) state at x = for the symmetric model (ao = 0) with 5 = 0. The 
light colored curves are the nuUclines / = and g = and the dark colored curves are 
the numerically computed trajectories. The trajectory through (0, 0) corresponds to the 
symmetric stationary front solution. The two other trajectories break the odd symme- 
try of the system and correspond to traveling fronts. The computational parameters are 
e = 1.0, (5 = 0, ai = 2.0, ao = for the stationary front and the same with e = .2 for the 
two traveling fronts. 

Figure 8: The front bifurcation line in the e — 6 plane for the symmetric (ao = 0) model 
with ai = 2.0. The solid circles represent the numerically calculated bifurcation line. The 
solid line is the linear approximation (valid for 5 <^ 1) and the dashed line is the singular 
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pertubation result (valid for 1). 

Figure 9: Phase portraits of front solutions connecting the {u+,v+) state at x = — 
to the (tt_,v_) state at x = for the symmetric model (ao = 0) with 5 > 1. The 
light colored curves are the nuUclines / = and ^ = and the dark colored curves are 
the numerically computed trajectories. The trajectory through (0, 0) corresponds to the 
symmetric stationary front with e > ec{S) and the two other trajectories correspond to the 
symmetry breaking traveling fronts with e < ec{S). 

Figure 10: (a) Relaxation of an initial pattern to the uniform down state for the regime 
of Ising fronts, e = 0.14, 5 = 0, ai = 2.0, ao = —0.1. (b) The convergence of the same ini- 
tial pattern toward a stable traveling pattern beyond the bifurcation from Ising to Bloch 
fronts, e = 0.09. 

Figure 11: Phase diagram in the e — ao plane. For the region to the right of the solid 
curve only one type of front solution exists and initial patterns do not persist. In the region 
between the solid and the dashed curves multiple stable fronts coexist but patterns still 
decay toward a uniform state. For the region to the left of the dashed curve initial patterns 
evolve toward persistent patterns in the form of stable traveling waves. Computational 
parameters are : S = 0, ai = 2.0. 

Figure 12: Time evolution of a front and a back (or two Bloch fronts following one 
another), (a) The back propagates faster than the front and binds to the front to form a 
traveling up state domain, (b) The back propagates slower than the front and the up state 
domain expands indefinitely. 

Figure 13: The onset of stationary domains in the Ising front regime as 5/e is increased. 
Ising fronts either (a) collide and anihilate to form a uniform state for e > fJ,stS, or (b) 
slow to a stop and form a stationary domain for e < fist^- 

Figure 14: Oscillating or breathing domains, (a) Steady oscillations close to the Hopf bi- 
furcation. Computational parameters e = 0.03,5 = 2.5, ai = 2.0, ao = —0.1. (b) Collapse 
of an oscillating domain further away from the Hopf bifurcation. Compuational parame- 
ters are the same with e = 0.025. 

Figure 15: The effect of the front bifurcation on the dynammics of two fronts propa- 
gating toward one another: (a) Below the bifurcation (e > e^S) an oscillating domain is 
formed. Computational parameters: e = 0.03,5 = 2.5, ai = 2.0, ao = —0.01. (b) Beyond 
the bifurcation the two fronts rebound from one another and propagate to the boundaries. 
Same parameters as in (a) except that e = 0.012. 

Figure 16: The effect of the front bifurcation on the dynamics of two fronts following one 
another: (a) Below the bifurcation (e > edS)) an oscillating domain is formed. Computa- 
tional parameters are e = 0.030, 5 = 2.5, ai = 2.0, ao = —0.012. (b) Beyond the bifurcation 
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a traveling domain is formed. Same parameters as in (a) except that e = 0.25. 
Figure 17: A typical bifurcation diagram for single domain structures. The solutions 
shown with the solid (dashed) line are stable (unstable) structures. The stationary domain 
solution loses stability to an oscillating domain at e = e;,,.. At e = etr, a branch of stable 
traveling domain solutions appears. By the symmetry c — > — c, x ~^ ~Xj this diagram is 
symmetric with respect to the c = axis, but only the positive speed branch of traveling 
domain solutions is shown. Also, there is an additional branch of zero speed solutions that 
is unstable for all values of e. Parameters: 5 = 2.5, ai = 2.0, ao = —.1, period = 100. 
Figure 18: Phase diagram for single domain structures in the e — 5 plane. In region I, 
domain structures are transient and the asymptotic state is either the uniform up or down 
state. In region II stable stationary domains coexist with the uniform up and down states. 
In region III domains oscillate with a typical region of steady oscillations denoted by the 
hatched area. In region IV domains travel. The boundary between regions III and IV (i.e. 
the onset of traveling domains) coincides with the front bifurcation line denoted by the 
thick solid/dashed curve. For this phase diagram, ai = 2.0 and ao = —0.1. 
Figure 19: Coexistence of traveling and oscillating periodic patterns beyond the bifurca- 
tion from Ising to Bloch fronts. Both (a) and (b) are obtained for the same computational 
parameters with different initial conditions. Note the difference in average wavelength 
between the two patterns. Computational parameters: e = 0.013, S = 2.5, ai = 2.0, ao = 
-0.1. 

Figure 20: Front-boundary interactions, (a) Transient domain region: the front is ab- 
sorbed at the boundary, (b) Stationary domain region: the front is stopped at the bound- 
ary, (c) Oscillating domain region: the front oscillates near the boundary, (d) Traveling 
domain region (beyond the front bifurcation): the front is reflected at the boundaries. All 
simulations are with Neumann boundary conditions. 
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Figure 1: Three basic cases for the nuUchnes of the model system (1.1). Case (i): The 
nuUchnes intersect at a single point on an outer branch of the cubic nuUcline. Case (ii): The 
nuUclines intersect at a single point on the middle branch of the cubic nuUcline. Case (iii): 
The nuUclines intersect at three points each on a different branch of the cubic nuUcline. 
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Figure 2: Patterns near and far from onset show the development of two spatial scales, (a) 
Stationary periodic pattern near onset of the Turing instability, (b) Stationary periodic pat- 
tern far from onset of the Turing instability, (c) Traveling pattern near the Hopf bifurcation, 
(d) Traveling pattern far from the Hopf bifurcation. In all the figures the solid line is the u 
field and the dashed line is the v field. 
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Figure 3: The free energy density S (see (3.3)) for different values of v. 
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Figure 4: Multiplicity of stable front solutions when e ^ 1. (a) Preparing the system in 
the down state and perturbing it locally at the left edge of the domain induces a front 
propagating to the right, (b) Preparing the system in the up state and perturbing it locally 
at the right edge of the domain induces a front propagating to the left. 
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Figure 5: Bifurcation diagrams of front solutions. The dots are data points representing the 
speed of the different types of stable front solutions that exist for each value of e. (a) The 
symmetric case (ao = 0). The solid line is the theoretical bifurcation diagram computed 
from (3.2f). (b) The nonsymmetric case (ao = .f) shows the unfolding of the pitchfork to a 
saddle node bifurcation. 
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Figure 6: The propogation direction of Bloch fronts is determined by the translation of the 
V field relative to the u field; v always lags behind u unless the front is stationary, (a) Left 
traveling front, (b) Stationary front, (c) Right travling front. 
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Fi gure 7: Phase portraits of front solutions connecting the state at 

X = -oo to the 

state at X = oo for the symmetric model (ao = 0) with (5 = 0. The light colored 
curves are the nuUclines / = and = and the dark colored curves are the numerically com- 
puted trajectories. The trajectory through (0, 0) corresponds to the symmetric stationary 
front solution. The two other trajectories break the odd symmetry of the system and corre- 
spond to traveling fronts. The computational parameters are e = 1.0, (5 = 0, ai = 2.0, ao = 
for the stationary front and the same with e = .2 for the two traveling fronts. 
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Figure 8: The front bifurcation line in the e — 6 plane for the symmetric (ao = 0) model 
with ai = 2.0. The solid circles represent the numerically calculated bifurcation line. The 
solid line is the linear approximation (valid for 6 <^ 1) and the dashed line is the singular 
pertubation result (valid for 6/e ^ 1). 
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Fi gure 9: Phase portraits of front solutions connecting the state at 

X = -oo to the 

state at X = oo for the symmetric model (ao = 0) with 6 > 1. The light colored 
curves are the nuUclines / = and g = and the dark colored curves are the numerically 
computed trajectories. The trajectory through (0, 0) corresponds to the symmetric station- 
ary front with e > €c{S) and the two other trajectories correspond to the symmetry breaking 
traveling fronts with e < €c{S). 
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Figure 10: (a) Relaxation ol an initial pattern to the uniform down state for the regime of 
Ising fronts, e = 0.14,(5 = 0,ai = 2.0, ao = —0.1. (b) The convergence of the same initial 
pattern toward a stable traveling pattern beyond the bifurcation from Ising to Bloch fronts, 
e = 0.09. 
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Figure 11: Phase diagram in the e — ao plane. For the region to the right of the soHd curve 
only one type of front solution exists and initial patterns do not persist. In the region between 
the solid and the dashed curves multiple stable fronts coexist but patterns still decay toward 
a uniform state. For the region to the left of the dashed curve initial patterns evolve toward 
persistent patterns in the form of stable traveling waves. Computational parameters are : 
,5 = 0, ai = 2.0. 
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Figure 12: Time evolution ol a Iront and a back (or two Bloch fronts following one another), 
(a) The back propagates faster than the front and binds to the front to form a traveling 
up state domain, (b) The back propagates slower than the front and the up state domain 
expands indefinitely. 
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Figure 13: The onset of stationary domains in the Ising front regime as 6/e is increased. 
Ising fronts either (a) coUide and anihilate to form a uniform state for e > fist^^ or (b) slow 
to a stop and form a stationary domain for e < fist^- 
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Figure 14: Oscillating or breathing domains, (a) Steady oscillations close to the Hopf bifur- 
cation. Computational parameters e = 0.03,(5 = 2.5, ai = 2.0, ao = —0.1. (b) Collapse of an 
oscillating domain further away from the Hopf bifurcation. Compuational parameters are 
the same with e = 0.025. 
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Figure 15: The effect of tlie front bifurcation on tlie dynammics of two fronts propagating 
toward one another: (a) Below the bifurcation (e > €c{S) an oscillating domain is formed. 
Computational parameters: e = 0.03,(5 = 2.5, ai = 2.0, ao = —0.01. (b) Beyond the bifur- 
cation the two fronts rebound from one another and propagate to the boundaries. Same 
parameters as in (a) except that e = 0.012. 
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Figure 16: The effect of tlie front bifurcation on tlie dynamics of two fronts following one 
another: (a) Below the bifurcation (e > €c{S)) an oscillating domain is formed. Computa- 
tional parameters are e = 0.030, 6 = 2.5, ai = 2.0, ao = —0.012. (b) Beyond the bifurcation 
a traveling domain is formed. Same parameters as in (a) except that e = 0.25. 
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Figure 17: A typical bilurcation diagram lor single domain structures. The solutions shown 
with the solid (dashed) line are stable (unstable) structures. The stationary domain solution 
loses stability to an oscillating domain at e = e^r- At e = ej^, a branch of stable traveling 
domain solutions appears. By the symmetry c — )• — c, x ~^ this diagram is symmetric 
with respect to the c = axis, but only the positive speed branch of traveling domain 
solutions is shown. Also, there is an additional branch of zero speed solutions that is unstable 
for all values of e. Parameters: 6 = 2.5, ai = 2.0, ao = —.1, period = 100. 
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Figure 18: Phase diagram for single domain structures in the e — 6 plane. In region I, 
domain structures are transient and the asymptotic state is either the uniform up or down 
state. In region II stable stationary domains coexist with the uniform up and down states. 
In region III domains oscillate with a typical region of steady oscillations denoted by the 
hatched area. In region IV domains travel. The boundary between regions III and IV (i.e. 
the onset of traveling domains) coincides with the front bifurcation line denoted by the thick 
solid/dashed curve. For this phase diagram, ai = 2.0 and ao = —0.1. 
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Figure 19: Coexistence of traveling and oscillating periodic patterns beyond the bifurcation 
from Ising to Bloch fronts. Both (a) and (b) are obtained for the same computational pa- 
rameters with different initial conditions. Note the difference in average wavelength between 
the two patterns. Computational parameters: e = 0.013,(5 = 2.5, ai = 2.0, ao = —0.1. 
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Figure 20: Front-boundary interactions, (a) Transient domain region: the front is absorbed 
at the boundary, (b) Stationary domain region: the front is stopped at the boundary, 
(c) Oscillating domain region: the front oscillates near the boundary, (d) Traveling domain 
region (beyond the front bifurcation): the front is reflected at the boundaries. All simulations 
are with Neumann boundary conditions. 



